A Spectral Approach to the Relativistic Inverse Stellar Structure Problem 
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A new method for solving the relativistic inverse stellar structure problem is presented. This 
method determines a spectral representation of the unknown high density portion of the stellar 
equation of state from a knowledge of the total masses M and radii R of the stars. Spectral 
representations of the equation of state are very efficient, generally requiring only a few spectral 
parameters to achieve good accuracy. This new method is able, therefore, to determine the high 
density equation of state quite accurately from only a few accurately measured [M, R] data points. 
This method is tested here by determining the equations of state from mock [M, R] data computed 
from tabulated "realistic" neutron-star equations of state. The spectral equations of state obtained 
from these mock data are shown to agree on average with the originals to within a few percent 
(over the entire high density range of the neutron-star interior) using only two \M, R] data points. 
Higher accuracies are achieved when more data are used. The accuracies of the equations of state 
determined in these examples are shown to be nearly optimal, in the sense that their errors are 
comparable to the errors of the best-fit spectral representations of these realistic equations of state. 

PACS numbers: 04.40.Dg, 97.60.Jd, 26.60.Kp, 26.60.Dd 
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I. INTRODUCTION 

The standard stellar structure problem consists of de- 
termining the structure of a star by solving the coupled 
gravitational and hydrodynamic equations with an as- 
sumed equation of state for the stellar matter. The so- 
lutions to the standard problem determine the various 
observable properties of the stars with a given equation 
of state like their total masss M, their total radii R, etc. 
The inverse stellar structure problem determines what 
equation of state is required to produce stellar models 
having a given set of macroscopic obserables. The goal of 
this paper is to find efficient and robust methods of solv- 
ing this inverse stellar structure problem. The method 
developed here is based on the use of spectral expan- 
sions to represent the equation of state. The values of 
the spectral coefficients in these expansions are fixed in 
this method by matching stellar models based on these 
equations of state to observed values of the masses and 
radii of the stars. Once fixed, these coefficients determine 
the equation of state that represents the (approximate) 
solution to the inverse stellar structure problem. 

For non-rotating stars in general relativity theory, the 
simplest version of the stellar structure equations were 
first derived by Oppenheimer and Volkoff 
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where m(r) represents the mass contained within a 
sphere of radius r; p(r) is the pressure; and e(p) is the 
total energy density of the fluid. Solving these equations 
with a given equation of state is the standard relativis- 
tic stellar structure problem. Once an equation of state 
e = e(p) is specified, these equations determine a one 
parameter family of stellar models, m = m(r,p c ) and 



p = p(r,p c ), where p c is the value of the pressure at the 
center of the star r = 0. These solutions then deter- 
mine various macroscopic properties of the stars, includ- 
ing their outer radii R{p c ) where p[R(p c ),p c ] = 0, and 
their total masses M(p c ) — m[R(p c ),p c ]. These macro- 
scopic properties are (at least in principle) observable. 

The standard stellar structure problem can be thought 
of as a map that takes the equation of state [a curve in 
energy density - pressure space e = e(p)], into a curve in 
the space of macroscopic observables, e.g. [M(p c ), R(p c )]. 
The inverse stellar structure problem consists of finding 
the inverse of this map Q , he. determining the equation 
of state of the stellar matter from a knowledge of some 
information about the macroscopic structures of the stars 
(like their masses M and radii R). The solution to this 
problem, like the solutions to many inverse problems, is 
less straightforward than the solution to the standard 
problem. 

The inverse stellar structure problem is probably more 
relevant for practical relativistic astrophysics, however, 
than the standard problem. The highest density part 
of the equation of state in neutron stars, for example, 
is not well known. Matter in this state is well beyond 
the reach of laboratory experiments, so there is no in- 
dependent way of directly measuring its properties, in- 
cluding its equation of state. Numerous attempts have 
been made to model this matter theoretically, but even 
today there is no consensus among theoreticians. Pre- 
dictions of the energy density for a typical neutron-star 
central pressure, for example, often differ by an order 
of magnitude. Therefore, the standard stellar structure 
problem for neutron stars is not terribly useful. In con- 
trast, the inverse problem may provide an important tool 
for learning about high density nuclear matter. Numer- 
ous high quality measurements of the masses of neutron 
stars are now available @, and a few (fairly imprecise 
and model-dependent) radius measurements are starting 
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to become available as well |j-|7[. In principle then, the 
inverse stellar structure problem should (eventually) al- 
low us to measure the high density equation of state of 
neutron-star matter directly. This measurement would 
provide important information about nuclear interactions 
that can not be obtained at present in any other way. 

One naive approach to solving the inverse stellar struc- 
ture problem for neutron stars would simply be to match 
their observed properties, e.g. their [M, R] data, with 
models of those stars based on different micro-physical 
models of the dense material in their cores. In this ap- 
proach the model equation of state whose stellar models 
best matches the data would be declared the observed 
neutron-star equation of state. This approach would 
clearly be ideal if there were wide consensus on exactly 
what the high density core material is, and if there were 
a reasonably simple model for this material that were 
known, up to a few undetermined parameters that could 
be fixed by these observations. Unfortunately the wide 
diversity of "realistic" neutron-star equations of state in 
the literature, suggest that (in the near term at least) 
this approach is not likely to be effective or conclusive. 

A more practical variation of this approach uses some 
knowledge about the expected properties of neutron-star 
matter in an intermediate range of densities, and a more 
empirical description of the equation of state for larger 
densities [H Q . This approach is more promising, but the 
proposed model equations of state of this type have many 
free parameters that must all be fit by the observational 
data. Since these data are likely to be sparse for some 
time, wc take a different approach here. 

Our goal is to find efficient and robust methods for 
solving the inverse stellar structure problem that use no 
prior knowledge of the high density micro-physics at all. 
A (somewhat impractical) method for solving the inverse 
stellar structure problem that uses no information about 
the micro-physics of the high density equation of state 
was given in the literature about 20 years ago @]. This 
traditional method can be summarized as follows. The 
total masses M and radii R of all of the stars associ- 
ated with a particular equation of state are assumed to 
be known. The equation of state is also assumed to be 
known up to some pressure pi with corresponding energy 
density = e(pi). Let Mi = M(jpi) and Ri = R(pi) de- 
note the mass and radius of the star whose central pres- 
sure is p c = pi. Now choose another point, [Mj+i,i?j+i], 
along the mass-radius curve, with a slightly larger central 
pressure. The outer layers of this new star are composed 
of low pressure material, p < Pi, where the equation of 
state is known. The stellar structure equations, (TTJ) and 
([2]) , can therefore be solved in this outer region starting at 
the surface of the star, r = Ri+i, where p(i?^+i) = and 
m(i?i+i) = A/i+i, by integrating inward toward r = 0. 
This integration can be continued until the point r = r; + i 
where p(r,+i) = pi and the known equation of state ends. 
This integration determines the radius r, +1 and the mass 
m i+i = m { r i+i) of a "small" core of high pressure mate- 
rial, p > pi where the equation of state is not yet known. 



If this core is small enough, the stellar structure equations 
can be solved there as a power series expansion about 
r = 0. The coefficients in this expansion are functions 
of the central density ej+i and pressure Pi+i of this little 
core. Since the mass and radius of this core are known, 
rrij+x and rj+i, this power series can be "inverted" to 
determine e%-\-\ and Pi+i 0. This new point [^4.1,^+1] 
provides a small extension of the equation of state beyond 
[ei,pi\. Iterating these steps then determines a sequence 
of closely-spaced points along the high density portion of 
the equation of state curve. 

This traditional solution to the inverse stellar structure 
problem is unfortunately very impractical. A large num- 
ber of points [Mj,i?j] are needed from the mass-radius 
curve to achieve modest accuracy in the calculation of the 
corresponding points [ei,Pi] along the equation of state 
curve. Since [Mi,Ri] points are very difficult to mea- 
sure (at least for neutron stars) it is unlikely that there 
will ever be enough data to use this traditional method to 
determine the unknown high density part of the neutron- 
star equation of state. 

This paper proposes a rather different approach to the 
solution of the inverse stellar structure problem, an ap- 
proach that can be very effective even when only a small 
number of [Mi,Ri] data points are available. The equa- 
tion of state in this new approach is expressed as a para- 
metric equation, e.g. e = e(p,7fc), instead of a table of 
values [e,-,Pi] . The parameters jk are adjusted to give the 
best-fit approximation to a particular equation of state 
model. Parametric representations of this sort, based on 
spectral expansions, have been shown to be extremely 
efficient at representing the high density portions of "re- 
alistic" neutron-star equations of state [9j. Only a few 
non- vanishing are generally needed to achieve 1% ac- 
curacy in most cases. 

The basic idea of this new method for solving the in- 
verse stellar structure problem is to choose the equation 
of state parameters 7& by minimizing the differences be- 
tween the masses and radii of real neutron stars, Mj 
and Ri , with those based on the parametric model equa- 
tion of state, M(p c ,7fc) and R(p c >lk)- Once the 7^ are 
fixed, the parametric equation e = e(p, 7^) then provides 
an approximate solution of the inverse stellar structure 
problem. Spectral expansions typically converge expo- 
nentially as the number of terms in the expansion are 
increased (for smooth functions). These approximate so- 
lutions to the inverse stellar structure problem are there- 
fore expected to converge to the exact equation of state 
as the number of data points [Mi, Ri] and the number of 
parameters 7fc fixed by this method are increased. 

The remainder of this paper presents details on how to 
implement this new spectral approach to the solution of 
the inverse stellar structure problem, along with practical 
tests of its accuracy and efficiency. Section [TT] reviews the 
particular spectral representation of the equation of state 
used in the solution presented here. Section UTTl describes 
how the spectral parameters 7^ are fixed by matching to 
the given [Mi,Ri] data points. Section HVl presents a se- 
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ries of numerical tests of the accuracy and efficiency of 
this new method. Mock [M,,i?j] data computed from a 
collection of 34 "realistic" neutron-star equations of state 
are used as input in these tests. These tests show, for ex- 
ample, that the resulting spectral equation of state agrees 
with the exact to within a few percent (on average) when 
only two [Mi, Ri] data points are used. Higher accuracies 
are (generally) achieved when more data points are used. 
Section [V] discusses some of the limitations of the numer- 
ical tests presented here, and proposes several ways that 
the basic method developed here might be extended and 
improved. Some of the more complicated technical de- 
tails needed to implement this method are described in 
two Appendices. Appendix lAl describes how to evaluate 
the derivatives of M(h c ,^k) and R(h Cl ^k), with respect 
to the parameters h c and 7^ . Appendix [B] describes the 
interpolation method used here to bridge the gaps be- 
tween points in the exact "realistic" equation of state 
tables. 



II. SPECTRAL REPRESENTATIONS OF THE 
EQUATION OF STATE 

The version of the relativistic stellar structure equa- 
tions most useful for our analysis here requires that the 
equation of state be written in a form where the energy 
density e and pressure p are given as functions of the rel- 
ativistic enthalpy, h. The usual form of the equation of 
state, e = e(p), must therefore be re- written as a pair of 
equations e = e(h) and p = p(h), where the enthalpy h is 
defined as 
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The needed expressions, e = e(h) and p — p(h), can 
be constructed by inverting h = h(p) from Eq. ([3]) to 
obtain p = p(h), and then by composing the result with 
the standard form of the equation of state, e = e(p), to 
obtain e(h) = e[p(h)]. 

The transformations needed to construct e = e(h) and 
p = p(h) in this way are difficult to perform efficiently 
and accurately in numerical computations. Therefore it 
is best to construct a spectral representation of the equa- 
tion of state that is based directly on h. This can be done 
by introducing an enthalpy based spectral expansion of 
the adiabatic index r @: 
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where ho is the lower bound on the enthalpy, ho < h, in 
the domain where the spectral expansion is to be used. 
This is a standard spectral expansion of the function 
logT(h) in which the [log(h/ho)] h are the spectral basis 



functions and the 7fc are the spectral expansion coeffi- 
cients (or parameters). 

The functions p(h) and e(h) are obtained from T(h) by 
integrating the system of ordinary differential equations, 
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that follow from the definitions of h and T in Eqs. ([3]) 
and ([4]). The general solution to these equations can be 
reduced to quadratures: 
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The constants po and eo are defined by po = 2?(^o) an d 
e = e(h ) respectively. While these quadratures can 
not be done analytically for the spectral expansion of 
r(/i) given in Eq. ([5]), they can be done numerically very 
efficiently and accurately using Gaussian quadratures. 1 

The method of solving the inverse stellar structure 
problem proposed in Sec. HIT] is based on this spectral 
representation of the equation of state: e = e(h,jk) and 
p = p(h,-fk)- Any equation of state can be represented 
approximately in this way by using a finite number of 
spectral parameters 7^ in the expression, Eq. ([5]), for 
T(h). In analogy with other spectral expansions, the ac- 
curacy of these approximations are expected to converge 
exponentially (for smooth equations of state) as the num- 
ber of spectral coefficients is increased. Numerical tests 
that fit "realistic" neutron-star equations of state using 
this method Q are consistent with this expectation about 
the convergence of these expansions. 



III. FIXING THE SPECTRAL PARAMETERS 

The spectral approach to the inverse stellar structure 
problem fixes the equation of state parameters 7^ by 
choosing the values whose stellar models best match a 



1 The numerical accuracy (and hence the efficiency) of the numeri- 
cal integrations in Eqs. (JSJl and (1 1 Oi l can be improved significantly 
by changing integration variables from h to x = log(/i//io) before 
performing standard Gaussian quadratures. Our tests achieved 
accuracies about fO^ 11 for e(h, 7^.) and p(/i,7fc) with 10 Gaus- 
sian integration points using the x variable, compared to about 
10 — 3 accuracies for the same tests using the h variable. 
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collection of points [Mi,Ri] from the exact mass-radius 
curve. The process of fixing the 7^ in this way can be 
made more efficient numerically by using a somewhat 
non-standard version of the stellar structure equations. 
Therefore, we digress briefly here to review this alternate 
formulation. 

The Oppenhcimer-Volkoff version of the stellar struc- 
ture problem, Eqs. ([I) and ©, determines m and p as 
functions of r. That traditional approach has two incon- 
venient features: First the integration domain, [0, R], is 
only known after the fact when the surface of the star at 
r = R is found numerically by solving p(R) = 0. Second, 
the equation p(R) = that defines the surface of the star 
is somewhat difficult to solve numerically because dp/dr 
typically vanishes at r = R. These inconveniences can be 
avoided by transforming the equations into a form where 
m and r are determined as functions of the relativistic 
enthalpy h (see Ref. 0]): 
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Solving the equations in this form begins by specifying 
"boundary" conditions, m(h c ) = r(h c ) = 0, at the center 
of the star where h = h c , and then integrating toward the 
surface of the star where h = 0. The derivative dr/dh 
in Eq. (|12|) is non-zero and bounded at the surface of 
the star, so this formulation completely eliminates the 
problems associated with solving p(R) = to locate the 
star's surface. This version of the problem is also eas- 
ier to implement numerically because it is carried out on 
the domain [h c , 0], which is fixed before the integration 
is performed. The total mass and radius of the stellar 
model are obtained in this formulation simply by evalu- 
ating the solutions m(h) and r(h) at the surface of the 
star where h = 0: M = m(0) and R = r(0). More details 
about how to implement this formulation of the stellar 
structure problem are given in Ref. Q and in Appendix El 
of this paper. 

This alternate formulation of the stellar structure 
problem, Eqs. (|TT|) and (|T2"j) . requires that the equation 
of state be expressed in terms of the enthalpy, i.e. that 
e = e(h) and p = p(h) be provided. The spectral rep- 
resentations, e = e(h, 7*.) and p = p(h, 7/.), described 
in Sec. [TT] are therefore ideal for this purpose. The gen- 
eral solution to this form of the stellar structure prob- 
lem is a pair of functions of the form m(h, h c ,jk) and 
r(h, h c ,"fk)- These solutions arc determined uniquely by 
the parameter h c , the central enthalpy of the star, and 
7/c, the spectral parameters that determine the equation 
of state. The total mass M(h c ,jk) and radius i?(/i c ,7fc) 
associated with one of these stellar models are deter- 
mined from these solutions by M(h c ,"fk) = m(0, h c ,jk) 
and R{h Cl j k ) = r(0, h c ,j k ). 

The new method of solving the inverse stellar struc- 
ture problem, which we introduce here, fixes the 
values of the spectral parameters 7^ by minimizing 



the differences between the model mass-radius values 
[M(/i c ,7fc),i?(/i c ,7fc)] and points [Mj,i?j] from an exact 
mass-radius curve. Thus we fix the values of the jk by 
minimizing 
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with respect to each of the 7^.. Each of the N stars stellar 
models used in this fit has a central enthalpy, h l cl whose 
value is also needed (along with the 7^.) to construct 
the mass-radius values [AI(h l c ,jk),R(h l c ,jk)}- Since the 
h\ will not be known a priori, these additional parame- 
ters must also be determined as part of the fitting pro- 
cess. These parameters are fixed, therefore, by minimiz- 
ing X 2 (^c!7fc) w hh respect to variations in each of the 

In summary then, this new method of solving the in- 
verse stellar structure problem determines the equation 
of state by fixing the spectral parameters in a way that 
minimizes the differences between the model mass-radius 
values [M(hl, 7^.), R{h\, 7^)] and values [Mi,Ri] from an 
exact mass-radius curve. This minimization problem is 
equivalent to solving the -/V st ars equations 



dx 2 
dhi 
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and the N-y k (the number of spectral parameters) equa- 
tions 



(15) 



for the parameters h l c and 7^. Since the number of in- 
dependent Mi and Ri data values used in this fitting 
process is 2 Altars, it follows that the maximum number 
of spectral parameters that can be fixed in this way is 

A number of numerical methods for solving non-linear 
least squares problems such as Eqs. (|T4l) and (jl"5j) are dis- 
cussed in the literature. Many of these methods require 
only that X 2 (^c;7fc) b e provided numerically for arbi- 
trary values of the parameters h % c and jk ■ Some methods 
also require in addition that the values of the deriva- 
tives d\ 2 /dh l c and dx 2 /djk be provided. The numerical 
tests described in Sec. IIVI use the Levenberg-Marquardt 
method for solving these equations, and this method re- 
quires that both the values and the derivatives of x 2 be 
provided. 

The derivatives of % 2 are determined by the deriva- 
tives of M(/i*,7fe) and R(h l c ,jk) with respect to h\ and 
7fe. These derivatives can be approximated numerically 
by expressions of the form, dMj djk ~ [M(h l c , 7fc + ^7^) — 
M(/i*,7fe — 6jk)]/2S-fk- We find it is more efficient (and 
more accurate) however to evaluate these derivatives by 
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solving an auxiliary system of ordinary differential equa- 
tions, that are obtained by differentiating Eqs. ([lip and 
(Ti"2|) with respect to these parameters. This method 
of evaluating the needed derivatives of M(/i*,7fc) and 
R(h l c ,jk) is described in some detail in Appendix El 



IV. TESTING THE SPECTRAL INVERSION 
METHOD 

In this section we test the spectral method of solving 
the inverse stellar structure problem by analyzing sets of 
mock [Mi, Ri] data points from mass-radius curves based 
on known "realistic" neutron-star equations of state. We 
use these mock [Mi, Ri] data to construct best-fit values 
for the spectral equation of state parameters 7^, using 
the least-squares method outlined in Sec. IIIII Then we 
compare the equation of state e(h, 7^) constructed in this 
way with the exact equation of state e(h) used to compute 
the mock [Mi,Ri] data points. We perform these com- 
parisons using different numbers of [Mi, Ri] data points 
to determine how the accuracy of the approximate equa- 
tion of state improves as the number of data points is 
increased. We construct the mock [Mi,Ri] data points 
using 34 different realistic neutron-star equations of state 
to explore how well the method works for a fairly wide 
variety of different equations of state. 

The 34 equations of state used to construct the [Mi, Ri] 
data points used in these tests are the same ones used by 
Read, Lackey, Owen and Friedman [l(| in their study 
of approximate piecewise polytropic fits to the equation 
of state. These 34 realistic equations of state are based 
on a variety of different models for the composition of 
neutron-star matter, and a variety of different models for 
the interactions between the particle species present in 
the model material. Descriptions of these realistic equa- 
tion of state models, and references to the original pub- 
lications on each of these equations of state are given 
in Ref. [Ioj |. and are not repeated here. The individual 
equations of state are referred to here using the abbrevia- 
tions used in Ref. [Ioj|, e.g. PAL6, APR3, BGN1H1, etc. 
The list of these equations of state are given in the first 
column of Table III of Ref. [l(| as well as the first column 
of Table U here. Spectral fits have already been shown to 
provide excellent approximations to these 34 equations of 
state in Ref. Only two or three spectral coefficients 
7fc are needed to achieve accuracies at the few percent 
level. The new question being studied here, therefore, is 
whether the spectral parameters 7^ can be determined 
by fitting [M,, R4] data instead of fitting directly to the 
equation of state itself. 

These 34 exact realistic equations of state are provided 
to us as tables of energy density and pressure points 
[e»,Pi]. We compute the enthalpy values hi corresponding 
to the points in these tables, interpolate between these 
tabulated points whenever necessary, and construct com- 
plete enthalpy based equations of state [e(h),p(h)] from 
these tabulated data using the methods described in Ap- 



pendix [Bj Whenever we refer to one of these exact real- 
istic equations of state, we mean the one constructed by 
interpolating the tabulated equation of state data in this 
way. 

We construct sets of mock [Mi, Ri] data points by solv- 
ing the standard stellar structure problem using each of 
the exact realistic equations of state described above. 
Figure [T] illustrates these exact mass-radius data for three 
of the realistic equations of state: PAL6, APR3, and 
BGN1H1. We select subsets of these points for each of 
our tests of the spectral inversion method. We limit the 
points used for our tests to small numbers of models 
(since we do not anticipate that large numbers of ob- 
servations are likely to be available any time soon) that 
fall within the astrophysically relevant range of masses: 
1.2M Q < M < M max , where M max is the maximum mass 
star allowed for a particular equation of state. We choose 
models for our tests that are (approximately) evenly 
spaced in mass within this range: 
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for i = 1, A r st ars- The large dots on each curve in Fig.Q] 
illustrate the points in the data sets with N s t ars = 5 from 
three of these equations of state. 
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FIG. 1: Exact mass-radius curves for the three realistic 
neutron-star equations of state, PAL6, APR3, and BGN1H1. 
Large dots illustrate the data points for the spectral method 
tests that use N ata , la = 5 stellar models. 

We have used these spectral methods to find solutions 
to the inverse stellar structure problem with mass-radius 
data sets [Mi,Ri] having Altars = 2, 3, 4 and 5 stellar 
models. In each case, we have constructed approximate 
equations of state with N lk non-vanishing spectral pa- 
rameters, for N lk = 2, Altars- We use the Levenberg- 
Marquardt algorithm for solving the non-linear least 
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squares problem, Eqs. (| 14[) and (fl5j) . as described in 
Ref. [ll| . We find that this method is extremely efficient 
at finding solutions that minimize the X 2 (^ci7fe) defined 
in Eq. (|13p. given a reasonably accurate initial guess for 
the values of the parameters h\ and jk ■ Our purpose here 
is to explore the overall accuracy of the spectral inver- 
sion method, and is not (at this point) directly aimed at 
producing an optimal robust method for analyzing real 
neutron-star observational data. Therefore we use our 
knowledge of the exact equation of state to provide the 
initial guesses for h l c and 7fc in the test solutions reported 
here. In particular we use the values of h c for the stel- 
lar models computed from the exact equation of state as 
initial guesses for h l c . Similarly we use the values of 7fc ob- 
tained by directly fitting the exact equation of state data 
as the initial guesses for these quantities. While these 
initial guesses are not precisely the solutions to the non- 
linear least squares problem, Eqs. (fT4")l and (fT5)) , they are 
close enough that the Levenberg-Marquardt algorithm 
easily converges. 

We assess the accuracy of our solutions by evaluat- 
ing the differences between the approximate equation of 
state e(h, jk) produced by the inversion process, and the 
exact equation of state e(h) that was used to construct 
the [Mi,Ri] data. We measure these differences by con- 
structing the error measure, 
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The sum in Eq. (|17p is taken over the points, [e^, hi], from 
the tabulated realistic equation of state. Only the iV 0OS 
points that lie within the range of densities £q < ti < e m ax 
present in the neutron-star models associated with a par- 
ticular equation of state are used in this sum. TablcUlists 
the values of A jv 7 for each of the 34 realistic equations of 
state used in our tests. The results reported in Tableware 
for tests that fit the maximum number of spectral param- 
eters, N lk = Altars, for each set of [Mi,Ri] data points. 
Fitting Altars = 2 data points gives equations of state ap- 
proximations with average errors of only a few percent, 
A2(Average) = 0.040. Using larger numbers of [Mi,Ri] 
data points (generally) results in higher accuracy approx- 
imations to the equation of state, with average values of 
A 3 = 0.029, A 4 = 0.023 and A 5 = 0.017. Thus, the 
spectral approach to the inverse stellar structure prob- 
lem is capable of giving high accuracy measurements of 
the high density equation of state using only a very small 
number of [Mj,i2j] data points. 

Table Q] also contains two additional measures of the 
accuracy of our test solutions to the inverse stellar struc- 
ture problem. One of these, 



- A r 



AT 
Af os ' 



(18) 



provides another way to measure the error in the ap- 
proximate spectral equation of state. The quantity AfJ R 



in Eq. (|18[) refers to the error in the approximate equa- 
tion of state obtained by fitting [M i; Ri] data, as defined 
in Eq. {I?). The quantity Af 05 in Eq. (HHJ) refers to 
the error of the best possible TV-parameter spectral fit to 
this particular equation of state. The values of A^ os 
for the equations of state studied here were determined 
in Ref. [9(, and are given in Table II of that reference. 
The quantity therefore measures the accuracy of the 
approximate spectral equation of state obtained by solv- 
ing the inverse stellar structure problem, relative to the 
accuracy of the best possible approximate iV-parameter 
spectral equation of state. Table U shows that (almost) 
all of these Tn measures are of order unity: the ap- 
proximate equation of state obtained with this spectral 
inversion method is almost as accurate as the best pos- 
sible TV-parameter spectral fit to the equation of state. 
Table U also contains the values of the quantity x(^cj Ik), 
defined in Eq. (fT3")) , for each of the test solutions found 
here. These values of xi^lilk) are all much less than 
unity, which shows that the least squares method is do- 
ing a good job of minimizing the differences between the 
model values of [M{h l c , 7^), R(h l c ,jk)] and the exact data 
points [Mi, Ri]. 

In addition to the accuracy measures given in Table U 
we have made in depth studies of the errors associated 
with a few of these solutions to the inverse stellar struc- 
ture problem. We have selected for closer examination 
the equation of state whose iV sta rs = 2 solution has the 
smallest error, PAL6 with A2 = 0.0034, the equation of 
state whose error is the median of the cases studied in 
these tests, APR3 with A 2 = 0.0266, and the equation 
of state having the largest error, BGN1H1 with A2 = 
0.1352. Figure [5] shows the quantity log[e(/i, 7fc)/e(/i)], 
which measures the fractional difference between the best 
fit model equation of state e(/i,7fe) and the exact PAL6 
equation of state e(h). This figure shows that the errors 
in these solutions to the inverse stellar structure problem 
become smaller as the number of data points used in the 
fits with N lk = Altars becomes larger. This figure also 
shows that the model equations of state e(h, 7fc) do a very 
good job of approximating the actual equation of state 
e{h) over the entire range of densities that are present in 
the interiors of neutron stars. Figures [3] and U illustrate 
the analogous error measures for the equations of state 
obtained from [Mi,Ri] data based on the APR3 and the 
BGN1H1 equations of state. 

These three cases illustrate the range of errors obtained 
using the spectral inversion method: the best case, a typ- 
ical average case, and the worst case. We point out that 
the worst case, BGN1H1, equation of state has a strong 
phase transition within the neutron-star density range. 
Non-smooth equations of state of this type are difficult 
to fit using spectral methods, and convergence in these 
cases is typically a power law rather than exponential. 
Many more spectral parameters are therefore needed to 
approximate these cases accurately. We point out, how- 
ever, that the values of Ajv for the BGN1H1 case are 
about ten percent, so even in this worst case the spectral 
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TABLE I: Accuracies of the neutron-star equations of state obtained by solving the inverse stellar structure problem. Ajv 
measures the average fractional error of the equation of state obtained by fitting N different [Mi, Ri] data pairs. The parameter 
Tjv measures the ratio of Am to the accuracy of the optimal iV-parameter spectral fit to each equation of state. The parameter 
Xn measures the accuracy with which the model masses M and radii R produced by the best fit equation of state match the 
exact data Mi and Ri. 
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inversion method gives a reasonably accurate estimate 
of the equation of state. The values of the Tjv for the 
BGN1H1 case all have values below 3.5, which shows that 
while it is difficult to model an equation of state of this 
type using a spectral fit, the spectral inversion method 
nevertheless does provide a solution that is comparable 
to the optimal Nj k -parameter spectral fit. 

Given an approximate spectral equation of state, 
e(/i, 7fc), we can use it to compute the complete mass- 



radius curve [M(h c , 7fc), R(h c , 7fe)] for the full range of 
central enthalpies, h c . This model mass-radius curve 
should agree with the exact curve [M{h), R(h)} very well, 
at least at the points [Mi,Ri] used in the inversion pro- 
cess. However, they will not agree everywhere, and the 
size of the differences is another measure of how well 
the approximate equation of state agrees with the exact. 
Figure [5] illustrates the differences between the model 
masses M(h c ,7k) arid the exact masses M(h c ) for the 
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FIG. 2: Ratios between various approximate equations of 
state, e(/i, 7fc), obtained by fitting [M, R] data, and the exact 
PAL6 equation of state, e(h). Note that log[e(/i, 7fc)/e(/i)] ~ 
[e(/i, 7fc) — e(/i)]/e(/i) measures the fractional error. 



FIG. 4: Ratios between various approximate equations of 
state, e(/i,7fc), obtained by fitting [M,R] data, and the ex- 
act BGNfNf equation of state, e(h). 
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FIG. 3: Ratios between various approximate equations of 
state, e(/i,7fc), obtained by fitting [M,R] data, and the ex- 
act APR3 equation of state, e(h). 



PAL6 equation of state. 2 Figures [5] and [7] make simi- 



lar comparisons for the APR3 and the BGN1H1 cases. 
We note that the error measures, log[M(h c ,jk)/M(h c )], 
shown in Figs. [5H71 are somewhat smaller in size than 
the error measures, log[e(h, 7fe)/e(/i)], shown in Figs.HHU 
These error measures provide one more piece of evidence 
that the spectral solutions to the inverse stellar structure 
problem are quite accurate. 



0.004 



& 0.003 



0.002 



O 



0.001 







PAL6 



N = 


N = 


2 


\ 


stars 




N = 


N, = 


3 


\ 


stars 




N = 


N, = 


4 


\ 


stars 




N = 


N = 


5 


\ 


stars 






1.5 2 2.5 

log( h c / h Q ) 



2 Note that the masses in Figs.[5]-[7]arc compared between models 
having the same central enthalpy h c . The central enthalpy of 
the exact model with M = Mi need not be exactly equal to the 
central enthalpy of the best fit model, h\ . Therefore these curves 
need not have zeros at those points. 



FIG. 5: Ratios between the masses, M(h c , 7fe) computed from 
the various fits, and M(h c ) computed from the exact PAL6 
equation of state, for a range of values of the central enthalpy 
h c of those models. 
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FIG. 6: Ratios between the masses, M(h c , 7^.) computed from 
the various fits, and M(h c ) computed from the exact APR3 
equation of state, for a range of values of the central enthalpy 
h c of those models. 
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FIG. 7: Ratios between the masses, M(h c ,^fk) computed from 
the various fits, and M(h c ) computed from the exact BGN1H1 
equation of state, for a range of values of the central enthalpy 
h c of those models. 



V. DISCUSSION 

In summary, we have developed a new method for 
solving the relativistic inverse stellar structure problem 
based on the construction of a spectral expansion of the 
unknown high density part of the equation of state of 
the star. The results of our numerical tests of this new 
method, described in Sec. IIV1 are quite impressive. Us- 



ing only two [Mi,Ri] data points, this new method can 
determine the entire high density part of the neutron- 
star equation of state with errors (on average) of just a 
few percent. The addition of more data points (gener- 
ally) results in higher accuracy approximations. We also 
show that TV-parameter spectral approximations to the 
equation of state determined in this way are almost as 
accurate as the best possible iV-parameter spectral ap- 
proximations. This is quite remarkable. It shows that 
macroscopic mass-radius measurements are strongly cor- 
related to the properties of the equation of state, and such 
measurements should therefore allow us (eventually) to 
measure the high density part of the neutron-star equa- 
tion of state with great precision. 

A close inspection of the results from the various tests 
summarized in Table U reveals a number of anomalies 
that merit further study. For example, the error measure 
A 7v Tfc , defined in Eq. ([TT]) , is expected to decrease as the 
number of spectral parameters N Jk is increased, i.e. that 
A at > Am 4-i . This seems to be true for most of our 
tests, but there are also a number of exceptions in Table|IJ 
The equation of state FPS, for example, has A2 = 0.0048, 
A 3 = 0.0061, A 4 = 0.0096, and A 5 = 0.0048. What is 
going on? Such a sequence of errors would be consistent, 
for example, with the idea that this particular equation 
of state is not well represented by these low order spectral 
expansions, i.e. that these expansions in this case are not 
yet in the convergent regime. This does not seem to be 
the case however since the optimal spectral fits to the 
FPS equation of state do appear to be convergent with 
these same numbers of spectral parameters, cf. Table 
II of Ref. @- Another (more likely) explanation of the 
anomalous results found in Table U is that the minima of 
X 2 (/i c ,7fc) found by the Levenberg-Marquardt algorithm 
for these cases are just local minima and not the desired 
global minima. An interesting area for further research 
on this problem, therefore, will be to explore the use of 
more robust numerical methods for finding global minima 
of complicated non-linear functions like x 2 (^c,7fc)- 

Another interesting direction for future research on 
this problem will be to explore how robust this kind of 
solution to the inverse stellar structure problem will be 
when applied to more realistic [Mi,Ri] data sets. The 
data used here were idealized in two important ways. 
First, the mock [Mi,Ri] data used here were supplied 
with very high precision. Real astrophysical measure- 
ments of these quantities will have significant errors. How 
will measurement errors influence the accuracy of the 
equation of state that is constructed by these techniques? 
Second, the mock [Mj,i?j] data used here were chosen 
to cover uniformly the astrophysically relevant range of 
neutron-star masses. Real astrophysical measurements 
will not be distributed in such an orderly way. How will 
the accuracy of the implied equation of state be affected 
by different, presumably less ideal, data distributions? 

The version of the inverse stellar structure problem 
studied here is based on the use of mass M and radius 
R measurements to determine the high density part of 
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the equation of state. These are not the only macro- 
scopic properties of neutron stars that could potentially 
be measured. It is not too difficult to imagine for ex- 
ample that the moment of inertias or the tidal Love 
numbers might be more easily observable for some types 
of observations. Another interesting direction for future 
study will therefore be to explore the use of other mea- 
surement data, say the mass and Love number (which 
could be measured using gravitational wave observations 
of neutron-star mergers), as input for solving the inverse 
stellar structure problem using the spectral methods de- 
veloped here. 
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2ir(e c + 3p c ) 
4(e c + 3p c ) 



1 1/2 



47T 

?ri5 = Anrf 



^3fc _ (fc +Pc) 2 
ri 5p c T c 



3(e c +Pc ) 
5p c r c 



2 1 



(A5) 
,(A6) 
(A7) 
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Appendix A: Computing Derivatives of M and R. 

This Appendix describes how the derivatives of the to- 
tal masses M(h c ,jk) and radii R{h c ,^k) are computed 
with respect to the parameters h c and 7^. To begin, 
however, we present a little more detail on how the alter- 
nate form of the stellar structure Eqs. (fTTj) and (fT2"]) are 
solved numerically. These equations are, 



dm 
~dh 
dr 
dh 



M(m, r, e,p) 
TZ(m, r,p) = - 



47rr 3 e(r — 2m) 
m + 4-7rr 3 p 
r(r — 2m) 
m + Airr 3 p ' 



(Al) 
(A2) 



The quantities e c , p c and r c in these expressions are the 
energy density, pressure and the adiabatic index evalu- 
ated at the center of the star where h = h c , e c = e{h c ), 
p c = p(h c ), and T c = T(h c ). 

It will be useful for our least-squares minimization 
problem to know how the solutions to Eqs. (|Aip and 
(|A2|) change as the parameters h c and jk are varied. Let 
A denote any one of the parameters: A = {h Cl "/k}- We 
wish to derive equations for the derivatives of the solu- 
tions to these equations with respect to these parameters: 
dm/dX and dr/dX. It is straightforward to determine the 
needed auxiliary equations by differentiating, Eqs. (|Aip 
and (|A2|) with respect to A: 



where the quantities A4(m,r, e,p) and lZ(m,r,p) merely 
represent the expressions on the right sides. These equa- 
tions are solved numerically by specifying conditions, 
m(h c ) = r(h c ) = 0, at the center of the star where h = h c 
and then integrating out to the surface of the star where 
h = 0. Like the standard Oppcnheimer-Volkoff version of 
the problem, Eqs. JT]) and ((2j), the right sides of Eqs. (|A1[) 
and (|A2|) . i.e. the functions Ai(m,r,e,p) and lZ(m,r,p), 
have the ill behaved form 0/0 there. Consequently it 
is necessary to start any numerical integration of these 
equations slightly away from the singular point h = h c . 
The needed starting conditions can be obtained using a 
power series solution to the equations. The needed power 
series are given in Eqs. (7) and (8) of Ref. 0, and can 



-( 


'dm\ 


dM dm 


dM dr 


dM de 


dh \ 


,~dXJ 


dm dX 


dr dX 


+ de dX 






dM dp 










dp <9A' 






d 




dTZ dm 


dTZ dr 


dTZ dp 


dh 




dm dX 


dr dX 


dp dX' 



(A9) 
(A10) 



The various derivatives dM/dm, etc. arc determined 
directly from the stellar structure equations, Eqs. (|A1[) 
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and 



A = 7^ these derivatives can be written as 



dM 
dm 
dM 

dr 
dM 

dp 
dM 
~d~T 

dU 

dm 

an 

dr 
dTZ 
dp 



87rr 3 e - M 
m + Anr 3 p ' 

2 3pM + 2e(2r - 3m) 



A-Kr 3 M 
m + Anr 3 p ' 
4irr 3 (r - 2m) 

m + 4nr 3 p 
2r-TZ 



m + 47rr 3 p ' 
l2nr 2 pK + 2(r - m) 
m + 47rr 3 p 



4 7rr 3p 



(All) 
(A12) 
(A13) 
(A14) 
(A15) 
(A16) 
(A17) 



For the case when A = 7^, the derivatives de/dj k and 
dp/djk are determined from Eqs. (|8l)- (fT0|) . The needed 
expressions are given by: 



dfi{h) 
djk 

d P (h) 
djk. 

de(h) 
djk 









f 


log 


{ho). 


I h 







dti 



-p(h) 



T(h>) ' 
h dfi(h') e h 'dti 



1 ho dlk [jl{h')f' 
dp(h) e(h) _ dp(h) e h p(h) 
dlk p{h) ~ dj k l^h)} 2 



(A18) 
(A19) 
(A20) 



The integrals needed to determine these quantities can 
be performed numerically with good efficiency and accu- 
racy using Gaussian quadrature. The equation of state 
does not depend on the parameter h c , and so de/dh c = 
dp/dh c = 0. Consequently the equations that determine 
dm/dh c and dr/dh c in Eqs. (g9| and (fATol) are some- 
what simpler than those for dm/djk and dr/dj k . 

The functions dm/dX and dr/dX are determined by 
solving Eqs. (|A9|) and (|A10[) numerically. This can be 
done by integrating them from the center of the star 
where h = h c out to the surface of the star where h = 0. 
To do this we need to impose the appropriate boundary 
conditions for these functions at h = h c . The needed 
boundary conditions can be found by differentiating the 
power series solutions, Eqs. (|A3| and (|A4[) . with respect 
to the parameters A. The quantities ri, r 3 , m 3 and 7715, 
which appear in these power series solutions, depend 
on the central values of the thermodynamic quantities 
e c = e(h c ), p c = p(h c ), and T c = T(h c ), and through 
them the parameters A = {/i c ,7fc}- For the case where 



dr(h) dri de c dri dp c 
djk [ de c dj k dp c dj k 

'dr^dec_ dr^dpc_ dr 3 dT c 
_de c dj k dp c d-y k dT c dj k 

+<D{h c -hf/ 2 , 
dm(h) 



{h c -hf' 2 

(h c h) 3 ' 2 



(A21) 



dj k 



dm 3 de c dm 3 dp c 



- h\ 3 / 2 



(h e - h) 



de c dj k dp c dj k _ 
dm§ de c dm$ dp c dm*, dT c 



de c d-f k dp c dj k 
+0(h c ~h) 7 / 2 . 



dT c d lk 



{h c -hf 2 
(A22) 



The derivatives of 77, r 3l m 3 and TO5 with respect to the 
parameters e c , p c and r c which appear in Eqs. (|A21[) and 
(|A22[> are given by: 



dri 

~dTc 
dri 
dp c 
dr 3 
de c 

dr 3 
dp c 

dr 3 
dT c 
dm 3 
de c 
dm 3 
dp c 

dm* 



2(e c + 3 Pc ) : 
dri 



de c ' 

r 3 dri 
ri de c 



4(e c + 3p c ) 

^3 dr! 3rj 
n dp c 4(e c + 3p c ) 

3ri(e e +p e ) 2 
20 Pc {e c + Sp^T 2 / 



4r 3 
ri 

4rg 

Tl 



(A23) 

(A24) 
6(e c + 3p c ) 



5p c r c 



5p 2 c T c 



4tt 



1 + 3e c dri 
ri de c 



dp c 



(A25) 

(A26) 
(A27) 

(A28) 

(A29) 



de c 



dm*, 
dp c 



dm*, 



Airrl 



'•3 



2e c r 3 dri 



n de 
A-Kr{{e c +p c ) 



dr 3 

' de c 



2r x +3(e c +Pc) 1 p- 
de c 



47re c r 1 



2r 3 dri dr 3 
ri dp c dp c 



4Trrf{e c +p c ) 



= 4-7rr 3 



e e dr 3 je c +p c ) 2 
n dT c 5p c Tl 



3p c (e c +Pc) dn 
ri dp c 



(A30) 

(A31) 
(A32) 



The values of the derivatives dp c /dj k and de c /d-f k are 
obtained by evaluating Eqs. (|A19|) and (|A20[) at h = h c , 
while the derivative dT c / d^f k is given by 



dr. 

dj k 



(A33) 
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For the case where A = h c the expressions for the 
derivatives dr/dX and dm/dX have somewhat different 
forms because h c appears explicitly in expansions in 
Eqs. (|A3[) and (|A4[) . Differentiating these series with 
respect to h c , keeping only the two leading terms, gives 



■{h c -K)- l t 2 



dr(h) n 
dh c ~ ~2 

dri de c dr\ dp c 
de c dh c dp c dh c 



3^3 
2 



+0(h 
dm(h) 3m3 



dhc 



2 

drnz de 



(h, 



- h) 1 ' 2 
dm 3 dp c 



5?Ti5 



de c dh c 
0{h c -hf' 2 . 



dp c dh c 



(h c 



(A34) 

hf/ 2 
(A35) 



The derivatives of n, r 3 , m 3 and m 5 with respect to 
the parameters e c and p c which appear in Eqs. (|A34[) 
and (|A35[) are given as before by the expressions in 
Eqs. (|S23]l - (|X3T]) . Th e deri vative s de c / dh c and dp c /dh c 
which appear in Eqs. (|A34|) and (|A35|) are obtained by 
evaluating Eqs. ([6]) and ([7J at h = h c : 



dp c 
dh c 
dec 
dh c 



(e c + Pcf 
PcF(hc) 



(A36) 
(A37) 



Appendix B: Interpolating the Exact Equation of 
State 



We are often presented with an "exact" equation of 
state that is represented as a table of energy densities 
Ej and the corresponding pressures pi. For our purposes 
here we will convert these to an equation of state of the 
form e = e(h) and p = p(h) in the following way. We 
do this by assuming that the exact equation of state is 
obtained for values intermediate between those given in 
the table, < e < ej+i, by the interpolation formula: 



V_ 

Pi 

Cj+l 



Ci + l 



\pg{p l+ l/pi) 

log(e i+ i/ej) 



(Bl) 
(B2) 



For smaller values of the density, e < ei, we assume: 



P_ 

Pi 

ci = 



log(P2/pi) 
log(e 2 /ei) ' 



(B3) 
(B4) 



Given this prescription for interpolation, it is straight- 
forward to show that the values of the enthalpy 



h{p) 



p dp' 



lo <P')+P r 
are given at the table entry values h% = h(jpi), by 



hi 

hi+i 



ci 



ci 

hi - 



■log 



ei +Pi 



Cj+i j q ej(ej+i +Pi+i) 
(k+i - 1 [ e l +i(e l +k) 



(B5) 

(B6) 
(B7) 



The pressure is determined as a function of the en- 
thalpy, by performing the integral in Eq. (|B5|) to give 
h{p), and then inverting. It is slightly more convenient 
to perform this inversion to give e(/i), from which it is 
straightforward to determine p(h) through Eqs. (|B3[) and 
dHU): 



e(h) = ei 



for h < hi, and 



Pi 



cxp ft 

Cl 



l/(ci-l) 



(B8) 



c(ft) 



e 4 +P t 
e. <! exp 



for hi < h < ftj+i. 



c,:+i - 1 



Ci+l 



l/(c i+1 -l) 

(B9) 
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